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Abstract 

The decay of very heavy metastable relics of the Early Universe can produce ultra- 
high energy cosmic rays (UHECRs) in the halo of our own Galaxy. On distance 
scales of the order of the halo size, energy losses are negligible — no Greisen-Zatsepin- 
Kuzmin cutoff is expected. In this letter we show that, as a consequence of the 
hierarchical build up of the halo, this scenario predicts the existence of small scale 
anisotropies in the arrival directions of UHECRs. We also suggest some consequences 
of this scenario which will be testable with upcoming experiments, as Auger. 



Heavy particles (mx ~ 10 12 — 10 14 GeV) can be produced in the Early Universe 
in different ways [1-4] and their lifetime can be finite though very long com- 
pared to the present age of the universe. The decay of these particles results 
in the production of UHECRs, as widely discussed in the literature [2,10-12]. 
If the relics cluster in galactic halos, as is expected, this can explain the cos- 
mic ray observations above ~ 5 x 10 19 eV. These particles may represent an 
appreciable fraction of the cold dark matter in the universe [1-4]. Can this 
class of models explain the details of the present data, including the spectrum, 
composition and anisotropy of the arrival directions? We will briefly review 
the status of the first two points, already discussed in the relevant literature 
and we will instead concentrate on the last point, paying particular attention 
to anisotropies one expects on small scales. 

The decay of heavy relics results usually in the production of a quark- ant iquark 
pair which rapidly hadronizes, generating two jets with approximately 95% of 
the energy in pions, and ~ 5% in baryons. The decay of the pions results in the 
observed high energy particles, mainly in the form of gamma rays, and in the 
generation of ultra-high energy neutrinos. The spectrum of the gamma pho- 
tons is relatively flat (~ E" 1 ' 5 ) reflecting the behaviour of the fragmentation 
function for the quarks. Therefore two main signatures of this model are: i) a 
flat energy spectrum; ii) composition dominated by gamma rays rather than 
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by protons. Moreover, as in all top-down models, heavy elements are expected 
to be completely absent. Unfortunately present data on the composition is 
extremely poor and it is impossible to rule out or confirm the presence of 
gamma photons in the UHECR events. 

In [14,11] the issue of the anisotropy was first addressed. The anisotropy re- 
sults from the asymmetric position of the Earth in the Galaxy, so that the 
flux of UHECRs coming from the direction of the galactic center should be 
appreciably larger than the flux from the anticenter direction. In [8,9] this 
issue was considered more quantitatively, taking into account the exposure of 
the present experiments. All authors concur that the present data is consistent 
with the predictions of the relic model for practically all reasonable values of 
the model parameters. 

Recently an interesting pattern has arisen from the analysis of the events with 
energy larger than 4 x 10 19 eV: in [5] the sample with this energy cut comprises 
47 events, whose overall distribution in space does not show appreciable devi- 
ation from isotropy. However, 3 doublets and one triplet were identified within 
an angular scale of 2.5°, comparable with the angular resolution of the exper- 
iment. A complete analysis, including the whole set of UHECR events above 
4 x 10 19 eV from the existing experiments was performed in [6]. This extended 
sample comprises 92 events and shows 12 doublets and two triplets (each 
triplet is also counted as three doublets) within an angle of 3°. The chance 
probability of having more than this number of doublets was estimated to be 
~ 1.5%. Although it is probably too soon to rule out the possibility that these 
multiplets are just a random fluctuation, it is instructive to think about the 
possibility that their presence contains some physical information about the 
sources of UHECRs. Clearly the most straightforward possibility is that the 
multiplets correspond to some local overdensity in the spatial distribution of 
the sources. Most of the top-down models for UHECRs (e.g. strings, necklaces, 
vortons, etc.) cannot naturally explain the multiplets. 

In the following we will discuss how the multiplets can be interpreted in the 
context of the super-heavy dark matter (SHDM) model. 

The existence of dark matter in galactic halos is all but established, and the 
main points are well understood in terms of hierarchical structure formation 
scenarios. High resolution N-body simulations, e.g., [7], suggest the following: 
%) the density of dark matter particles in galaxy size halos is peaked in the 
center; the density cusp scales as ~ r~ 7 with 7 ~ 1 — 1.5 on distances r which 
are much smaller than a core radius, which is of the order of several kpc in size; 
ii) outside the core, the slope of the profile steepens; it scales as oc r~ 3 at large 
distances; in) the profile is not completely smooth; some of the dark matter 
is in small clumps. To model the first two of these findings we adopt a dark 
matter density distribution in the form suggested by numerical simulations 
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where r c is the core size and n° is a normalization parameter. These two pa- 
rameters can be set by requiring that the halo contains a given total mass 
(Mjj) and that the velocity dispersion at some distance from the center is 
known (in the case of the Galaxy, the velocity dispersion is ~ 200 km/s in the 
vicinity of our solar system.). Alternative fits to the simulated dark matter ha- 
los and a discussion of whether or not simulated halos appear to be consistent 
with observations are provided in [7]. 

In addition to the smooth dark matter distribution, represented by eq. (1), N- 
body simulations also show that there is a clumped component which contains 
~ 10 — 20% of the total mass. The presence of these clumps are a natural 
consequence of the way in which gravity assembles dense virialized halos such 
as our galaxy today from the initially smooth density fluctuation field which 
was present when the cosmic microwave background (CMB) decoupled from 
the baryons. Simulations suggest that most of the mass which makes up a 
galactic halo was assembled by merging smaller clumps together at about 
z ~ 3. Much of the mass initially in a small clump which falls onto and orbits 
within the larger halo after that time gets tidally stripped from it. The amount 
of mass which is lost from any given clump increases as the distance of closest 
approach to the galactic halo center decreases; the mass is stripped away, from 
the outside in, as the clump falls towards the center. We will call the size of 
a clump, after its outsides have been stripped away, the tidal radius of the 
clump. Dynamical friction makes the clumps gradually spiral in towards the 
halo center. 

This has three main consequences. Firstly, the range of clump masses in the 
halo at the present time is different from that which fell in. Secondly, halos 
of a fixed mass do not all have the same tidal radius. Thirdly, the spatial 
distribution of the clumps in the halo is not the same as that of the dark 
matter. We found that a good fit to the joint distribution in clump mass and 
position in the simulations of [7] is 



where n° d is a normalization constant, rf is the core of the clumps distribution, 
and a describes the relative numbers of massive to less massive clumps. The 
simulations suggest that a ~ 1.9 [7]. The constraints on the core size are 
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weaker — we will study the range where rf is between 3 and 30 percent of Ru- 
in [7], a halo with M# « 2 x 10 12 M contains about 500 clumps with mass 
larger than ~ 10 8 M & . This sets the normalization constant in eq. (2). 

Clumps in the parent NFW halo are truncated at their tidal radii. The tidal 
radius of a clump depends on the clump mass, the density profile within the 
clump, and on how closely to the halo center it may have been. We assume 
that clumps of all mass are isothermal spheres (even though they are not truly 
isothermal, [7] suggest this is reasonably accurate): p c i(r d ) oc l/r d , where r d 
is the radial coordinate measured from the center of the clump. The tidal 
radius of a clump (R c i) at a distance r from the center of the parent halo is 
determined by requiring that the density in the clump at distance R c i from its 
center equals the local density of the NFW halo at the distance r. This means 
that 



where x c = r c /R H , x = r/R H and Rh is the virial radius of the halo, of order 
300 kpc. The average overdensity within R H is about 200 [7]. 

As shown in [11,9], the total (energy integrated) flux of UHECRs per unit 
solid angle from a smooth distribution of dark matter particles in the halo is: 



where R is the distance from the detector, and r is the distance from the 
galactic center (so R and r are related by trigonometrical relations accounting 
for the off-center position of the Earth in the Galaxy). The upper limit, R ma x, 
depends on the line of sight. 

The existence of a clumped component changes the flux in eq. (4) only in 
that uh should be replaced with the total dark matter density, the sum of the 
smooth and the clumped components. It is intuitively obvious that clumped 
regions will give an excess of events from certain directions, as was first pointed 
out in [13]. 

To see how important the clumped contribution is, we used two different ways 
of simulating the observed number of events. The first approach consisted of 
calculating the flux per unit solid angle [eq. (4)] along different lines of sight 
directly, taking into account the smooth plus clumped contributions to the 
total density profile [eqs. 1 and 2]. Once a smooth flux map distribution had 
been obtained, the UHECR events were generated from this distribution. In 
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the second approach, the events were generated in two steps. First, a random 
subset of the dark matter distribution, which is supposed to represent the 
subset of particles which decayed, was generated. The second step was to 
draw particles from this distribution, and then weight by the probability that 
the event would actually have been detected — so a chosen particle generates an 
event with probability oc 1/r 2 , where r is the distance between the particle and 
the detector. In both codes, the detector was assumed to be at the position of 
the Earth in the Galaxy, and a cut on the directions of arrival was introduced 
to account for the exposure of the AGASA experiment (taken here as an 
example). 

Fig. 1 shows as an example one of the generated flux maps: the map represents 
the ratio of the total flux including the contribution from clumps, to the flux 
obtained by using a smooth NFW profile. The various free parameters were 
r c = 8 kpc, rf = 10 kpc, and the mass distribution was truncated at a clump 
mass of 1% of the mass of the NFW halo. This sort of plot emphasizes the 
clump contribution. (To avoid confusion many of the smallest scale fluctua- 
tions, which arise from the many low mass clumps, have not been shown.) 

To calculate the small scale anisotropics, we generated 10 4 mock samples, 
each of 92 observed events, and counted the number of doublets and triplets 
for angular scales of 3, 4, and 5 degrees. Our codes can also be used to check 
the corresponding numbers for the case of isotropic arrival directions (as in 
[6]). Two sets of values of the cores for the NFW and the clumped component 
were adopted, one in which r c = 8 kpc and rf = 10 kpc (case 1) and the other 
with r c = rf = 20 kpc (case 2). The observed numbers of doublets within 3, 4, 
and 5 degrees for an isotropic distribution of arrival directions are given in [6] 
and are 12, 14 and 20 respectively. The number of doublets that we obtain in 
case 1 are 8, 14, and 21 within 3, 4, and 5 degrees respectively. The probability 
that the number of doublets equals or exceeds that observed is 12%, 47% and 
57% respectively. This should be compared with the 1.5%, 13.4% and 15.9% 
quoted in [6] for an isotropic distribution of arrival directions. 

We repeated the same calculation for the case 2. The corresponding averages 
and probabilities of exceeding the observed number of doublets within 3, 4 
and 5 degree scales are 6.6, 12, and 18, and 4.5%, 29% and 36% respectively. 

In both cases 1 and 2, the number of doublets on angular scales of 4 and 5 
degrees is consistent with the observed values; presumably the discrepancy at 
3 degrees is random chance. 

We have also studied the occurence of triplets. There is some ambiguity as to 
how a triplet is best defined; we have chosen to define triplets as configurations 
in which all three pairs would have been classified as doublets. (This means, 
for example, that a co-linear configuration of two doublets is not necessarily a 
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triplet.) With this definition, the average number of triplets in case one is 0.5, 
1.5 and 3, with the probability of having more than the observed triplets (2, 
2, 3 respectively) equal to 4%, 16% and 35%. For case 2, the correspondent 
numbers are 0.4, 1, and 2.5 triplets and 2%, 8% and 20% for the probabilities 
to have more triplets than observed. 

What is responsible for the multiplet-events in the SHDM model? If we study 
the case in which all the halo mass is in the smooth NFW component, then 
the number of doublets typically drops by one or two. This suggests that 
the anisotropy due to our position in an NFW halo can result in a number 
of multiplets of events which is considerably larger than if the arrivals were 
from an isotropic background. The number of multiplets from the clumped 
component is mainly affected by the presence of large nearby clumps, whose 
number depends on the high mass cutoff imposed in the mass function of 
clumps. A maximum mass of 1% of the halo mass implies a total mass in the 
clumps of ~ 10 — 15% of Mh, consistent with the results of the simulations 
[7]. Larger cutoffs imply larger mass fractions, which are harder to reconcile 
with the N-body simulations. 

Future experiments will definitely represent the real test for these sorts of 
models. Therefore it is useful to propose tests to be performed in the next 
years, in particular with the upcoming Auger experiment [17]. In the following 
we investigate what a full sky experiment would observe if UHECRs were 
produced by the decay of super-heavy halo dark matter. As discussed in [16] 
a powerful statistical tool that will be available with a full sky experiment is 
the angular power spectrum, as a function of I: 

cii) ~E «L (5) 

' m=—l 

where the coefficients a% m are defined as 

aim = -j^^2Yim(0,(f>), (6) 

and we assumed, for simplicity, that the exposure is flat over all the sky. Here 
N ev is the number of events and Y\ m are the real valued spherical harmonics, 
calculated in the direction determined by the two angles 6 and of the event. 

We simulated 5000 events over all the sky and computed the power spectrum 
in two cases: 1) purely NFW profile (no clumps); 2) full dark matter halo, 
(NFW plus a clumped component). The parameters were chosen as in case 
1 above. The calculation was done for one specific realization of clumps in 
the halo, since we are only interested here in the general features which may 
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appear in the power spectrum (the same realization shown in Fig. 1). 

The solid line (crosses) in Fig. 2 is the power spectrum corresponding to the 
smooth profile case (no clumps), while the dashed line (diamonds) is the power 
spectrum for the full dark matter profile (including clumps). In both spectra 
the dipole anisotropy is evident. This anisotropy is due to the asymmetric 
position of the Earth in the galaxy; it arises from the fact that the flux from 
the direction of the galactic center far exceeds (by an order of magnitude) that 
from the antigalactic center. As shown in [8,9] the present observations are not 
inconsistent with this predicted dipole anisotropy. Although not evident in this 
realization, there are cases in which the power spectrum (of the smooth plus 
clumped models) has some features (bumps). These are generally correlated 
with the presence of a nearby clump of dark matter which gives a substantial 
contribution to the flux on a specific angular scale. 

Contrary to the dipole enhancement, which is generic, the amplitude and 
frequency of the other bumps completely depends on the specific substructure 
distribution of the Milky Way's halo in which we happen to live. It is, therefore, 
difficult to quantify the likelihood of detecting these smaller bumps in the 
power spectrum of a future experiment. However, if such structures in the 
power spectrum are detected, then they are easily understood in the scenario 
discussed in this letter. 

In addition to features in the angular power spectrum (dipole and bumps) 
analyses of the composition of events will be precious in constraining or con- 
firming the scenario of SHDM in the halo; as discussed above, this should be 
dominated by gamma rays. For this reason, the halo magnetic field does not 
affect our conclusions. (If the proton and gamma ray fractions in UHECRs 
generated by SHDM are comparable [10], then the magnetic field might in- 
duce deflections on scales of a degree or so. Since we considered slightly larger 
scales, our conclusions are unlikely to change dramatically even in this case.) 

The study of the composition, together with an improved measure of the 
spectrum of UHECRs, should nail down the nature of the "real" sources of 
UHECRs and confirm or rule out the SHDM model. 

Aknowledgments This work was supported by the DOE and the NASA 
grant NAG 5-7092 at Fermilab. 



References 

[1] D.J.H. Chung, E.W. Kolb, and A. Riotto, Phys. Rev. Lett. 81, 4048 (1998). 



7 



[2] V.S. Berezinsky, M. Kachelriess, and A. Vilenkin, Phys. Rev. Lett. 79 4302 
(1997). 

[3] V.A. Kuzmin, presented at the Workshop "Beyond the Desert", Castle 
Ringberg, 1997, |astro-ph/9709187 ; International Workshop on Non- Accelerator 
New Physics, Dubna, 1997. 



[4] V.A. Kuzmin, and I.I. Tkachev, preprint |astro-ph/9903542 . 

[5] M. Takeda, Astrophys. J. 522, 225 (1999). 

[6] Y. Uchihori, et al., Astropart. Phys. 13, 151 (2000). 

[7] S. Ghigna, B. Moore, F. Governato, G. Lake, T. Quinn, J. Stadel, preprint 
|astro-ph/9910166t B. Moore, S. Ghigna, F. Governato, G. Lake, T. Quinn, J. 
Stadel, P. Tozzi, Astrophys. J. 524, 19 (1999); G. Tormen, A. Diaferio, D. Syer, 
MNRAS, 299, 728 (1998). 

[8] V.S. Berezinsky and A. Mikhailov, Phys. Lett. B 449, 237 (1999). 

[9] G.A. Medina Tanco and A. A. Watson, Astropart. Phys. 12, 25 (1999). 

[10] M. Birkel and S. Sarkar, Astropart. Phys. 9, 297 (1998). 

[11] V.S. Berezinsky, P. Blasi and A. Vilenkin, Phys. Rev. D58, 103515 (1998). 

[12] P. Blasi, Phys. Rev. D60, 023514 (1999). 

[13] V. Berezinsky, Invited talk at TAUP-99, Paris, September 6 - 10, 1999, preprint 
[hep-ph/0001163| . 

[14] S.L. Dubovsky and P.G. Tynyakov, Pis'ma Zh. Eksp. Teor. Fiz. 68, 99 (1998) 
[JETP Lett. 68, 107 (1998)]. 

[15] J.F. Navarro, C.S. Frenk and S.D.M. White, Astrophys. J. 462, 563 (1996). 



[16] P. Sommers, preprint astro-ph / 00040 1 6 



[17] J.W. Cronin, Rev. Mod. Physics 71, 175 (1999). 



8 



Figure Caption 



Fig. 1 Map of the ratio of the local flux of UHECRs obtained from a realistic 
dark matter distribution in the halo (NFW plus clumps) and the pure NFW 
case, for a specific realization of clumps in the halo. 

Fig. 2 Flux power spectrum for the NFW case (solid line with crosses) and the 
NFW plus clumped case (dashed line with diamonds). The same realization 
of clumps as in Fig. 1 has been used. 
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